% COZZI AND IMPULLITTI (2015), GLOBALIZATION AND WAGE POLARIZATION, REVIEW
% OF ECONOMICS AND STATISTICS

% CALIBRATION FUNCTION FILE 2: MOMENTS FUNCTION FILE

function Objective = Cali_fun_Piketty(X)


x0plus=[ 0.84542445670327   0.84266619012432   1];

global sig1 gam beta phi AD AF epsi fi1 n lambda eta rho eserv fi tau zf1 zd1


%======================== Externally calibrated ===========================


rho=0.03;      % discount factor

lambda=1.3;    % quality jump
n=0.0114;      % population growth rate
gam = 0.6;     % upper bound for share of skilled workers
eta = 1.5;     % Elasticity of substitution
fi1 = 0.0256;  % R&D difficulty parameter
V = 52;        % Life duration after college
TH = 4;        % Year of college
tau =1.74;     % Trade cost
zf1=1;         % Technology parameter production F
AF = 1;        % Technology parameter innovation F


% ================  Internally calibrated  ===============================


zd1 = X(1);    % Technology parameter production D
eserv = X(2);  % Service efficiency parameter
epsi = X(3);   % Ability distribution curvature
beta = X(4);   % Factor intensity production
phi = X(5);    % Factor intensity innovation
AD = X(6);     % Technology parameter innovation D



sig1=((exp(rho*V)-1)/(exp(rho*(V-TH))-1));
fi = ((exp(n*(V-TH))-1)/(exp(n*V)-1));

options=optimset('MaxFunEvals',100000000,'MaxIter',100000000,'TolFun',1e-7,'TolX',1e-7);


 [x,fval,exitflag,output]=fsolve(@CIcali_Piketty,x0plus,options);

 
 exitflag
% 
%             if exitflag<0
%                 disp('Problem here!');
% %                 pause
%             return
%             end


the0D = x(1);
the0F = x(2);
wDL = x(3);

wDH = (sig1*(the0D/(the0D-gam))*x(3));
wFL = (((the0F-gam)/(sig1*the0F))^(1-beta));
wFH = (sig1*(the0F/(the0F-gam))*wFL);

% ====================== TECHNOLOGY  ======================================

ZD = ((1/zd1)*(wDL^beta)*(wDH^(1-beta)));
ZF = ((1/zf1)*(wFL^beta)*(wFH^(1-beta)));
FD = ((1/AD)*((wDL^phi)*(wDH)^(1-phi)));
FF = ((1/AF)*((wFL^phi)*(wFH)^(1-phi)));

   % Marginal products: production
   
ZDL = ((1/zd1)*(beta*((wDH/wDL)^(1-beta))));
ZDH = ((1/zd1)*((1-beta)*((wDH/wDL)^(-beta))));
ZFL = ((1/zf1)*(beta*((wFH/wFL)^(1-beta))));
ZFH = (1/zf1)*((1-beta)*((wFH/wFL)^(-beta)));
   % Marginal products: R&D
   
FDL = ((phi*(1/AD)*((wDH/wDL)^(1-phi))));
FDH = (((1-phi)*(1/AD)*((wDH/wDL)^(-phi))));
FFL = ((phi*(1/AF)*((wFH/wFL)^(1-phi))));
FFH = (((1-phi)*(1/AF)*((wFH/wFL)^(-phi))));
% ========================== CUTOFFS ======================================

   % Skill cutoff to hire service workers
   
theHSD_fun = @(y) 1-y^epsi-(((wDH/wDL)*(1-eserv)*(y-gam))^epsi);
theHSF_fun = @(y) 1-y^epsi-(((wFH/wFL)*(1-eserv)*(y-gam))^epsi);

theHSD0 = 0.6;
theHSF0 = 0.6;

theHSD = csolve(theHSD_fun,theHSD0,[],1e-7,300);
theHSF = csolve(theHSF_fun,theHSF0,[],1e-7,300);

   % Skill cutoff for unskilled in production

theLSD = ((1-(theHSD^epsi))^(1/epsi));
theLSF = ((1-(theHSF^epsi))^(1/epsi));
  
   % Labor supplies
   
LDm = ((epsi/(1+epsi))*(the0D^(1+epsi)-theLSD^(1+epsi)));  % unskilled supply of production workers
LFm = ((epsi/(1+epsi))*(the0F^(1+epsi)-theLSF^(1+epsi)));

HD = ((1-the0D^epsi)*fi*(((epsi*((1-theHSD^(1+epsi))*(2 - eserv) + theHSD^(1+epsi) - the0D^(1+epsi)))/((1+epsi)*(1-the0D^epsi)))-gam));
HF = ((1-the0F^epsi)*fi*(((epsi*((1-theHSF^(1+epsi))*(2 - eserv) + theHSF^(1+epsi) - the0F^(1+epsi)))/((1+epsi)*(1-the0F^epsi)))-gam));

   % Average workers skills
   
theDm_avg = ((epsi/(1+epsi)) *((the0D)^(1+epsi)-(theLSD)^(1+epsi))/(the0D^epsi-theLSD^epsi)); %production
theFm_avg = ((epsi/(1+epsi)) *((the0F)^(1+epsi)-(theLSF)^(1+epsi))/(the0F^epsi-theLSF^epsi));


theDS_avg = ((epsi/(1+epsi)) *theLSD);  % Unskilled: Services
theFS_avg = ((epsi/(1+epsi)) *theLSF);

theDH_avg = (((epsi*((1-theHSD^(1+epsi)) + theHSD^(1+epsi) - the0D^(1+epsi)))/((1+epsi)*(1-the0D^epsi)))-gam);  % This average corrects for the increased supply of labour by skills: it is the average over the individuals of the averages per unit time, which means that those who work 2-eserv hours have their wages divided by (2-eserv). This is why we do not see (2-eserv) there, which of course we see it in the aggregate skill supply
theFH_avg = (((epsi*((1-theHSF^(1+epsi)) + theHSF^(1+epsi) - the0F^(1+epsi)))/((1+epsi)*(1-the0F^epsi)))-gam);  % This average corrects for the increased supply of labour by skills: it is the average over the individuals of the averages per unit time, which means that those who work 2-eserv hours have their wages divided by (2-eserv). This is why we do not see (2-eserv) there, which of course we see it in the aggregate skill supply


theDunskill_avg = ((epsi/(1+epsi))*the0D); % Average unskilled (service+production)

% ========== OTHER ENDOGENOUS VARIABLES ===============================

qF_til = ((1-fi1)*((lambda^(eta-1))-1)/n);
Fr = (FD/FF);
cD_til =( (Fr*(lambda-1)-(lambda-ZD*tau)) / ((lambda-ZD/tau)-(lambda-1)*Fr) );

x_til = ((1/FF)* (((lambda-1)/lambda)*(1+cD_til))/ (rho+(n/((1-fi1)*((lambda^(eta-1))-1))) + (n*fi1/(1-fi1))));

IF = (1/qF_til)*(1/(cD_til+1));
CF = ((1/IF)*(LFm / (qF_til*(ZFL/lambda)*(1+cD_til) + x_til*FFL)));

qF = (qF_til*IF);

qD = (1-(qF));

ID = ((1/qF_til)-IF);


CD = (qD/qF)*CF;

xindex = ((x_til)*CF); % R&D difficulty index




XX1 = CD;
XX2 = CF;
XX3 = ID;
XX4  = IF;
XX5  = wDL;
XX6  = wDH;
XX7  = wFL;
XX8 = wFH;
XX9  = xindex;
XX10  = wDH*theDH_avg; %avg skilled wages D
XX11  = wFH*theFH_avg; %avg skilled wages F
XX12  = wDL*theDm_avg;  %avg production wages D
XX13  = wFL*theFm_avg;  %avg production wages F
XX14 = wDL*theDS_avg; % avg service wages D
XX15  = wFL*theFS_avg; % avg service wages F
XX16  = ((wDH*theDH_avg)/(wDL*theDm_avg)); % skill premium D: skilled/production
XX16b = (wDH*theDH_avg)/(wDL*theDunskill_avg ); % skill premium D: skilled/unskilled
XX17  = ((wFH*theFH_avg)/(wFL*theFm_avg)); % skill premium F
XX18  = ((theDm_avg)/(theDS_avg)); % Production premium D
XX19  = ((theFm_avg)/(theFS_avg));  % Production premium F
XXstud = (1-the0D^epsi)*((exp(n*V)-exp(n*(V-TH)))/(exp(n*V)-1)); %share of students D
XX20  =(((1-(theHSD)^epsi)*fi)+((theHSD^epsi)-(the0D^epsi))*fi);  % share of skilled workers (physical units, not accounting for efficiency)
XX21 = ((((the0D)^epsi - (theLSD)^epsi))); % share of workers in production
XX22 = ((theLSD^epsi)); % share of workers in service
XX22bis = (((theLSD^epsi))+XXstud); % share of service sector workers plus students
GDPD = qD*(CD+CF) + (wDL*FDL + wDH*FDH)*xindex*ID; %wDL*(1-theHSD^epsi) + qD*(CD+CF) + (wDL*FDL + wDH*FDH)*xindex*ID; THE WAGES OF THE SERVICE WORKERS ARE PAID BY THE SKILLED WORKERS, SO THEY CANCEL OUT FROM AGGREGATE INCOME!  
GDPF = qF*(CD+CF) + (FFL + wFH*FFH)*xindex*IF; %wFL*(1-theHSF^epsi) + qF*(CD+CF) + (FFL + wFH*FFH)*xindex*IF; THE WAGES OF THE SERVICE WORKERS ARE PAID BY THE SKILLED WORKERS, SO THEY CANCEL OUT FROM AGGREGATE INCOME! 
XX24 = (ID/(ID+IF)); % Patent share
Total_population = XX20+ XX21+ XX22+XXstud;

x0plus=x;

CapitalD = FD*xindex*qD;% total value of the domestic firms - from eq (74)
CapitalF = FF*xindex*qF;% total value of the domestic firms - from eq (74)




betaPikD = CapitalD/GDPD; % Wealth-income ratio for country D
betaPikF = CapitalF/GDPF; % Wealth-income ratio for country F
 
% 
% Results = [XX6./XX5, XX8./XX7, XX5, XX7, the0D, the0F, theLSD, theLSF];
% disp('     wDH/wDl    wFH/wFL    wDL     wFL     the0D    the0F    theLSD     theLSF')
% disp(Results)

% ======================== CALIBRATED MOMENTS ============================

% Moments = [XX16(1,1), XX18(1,1), XX20(1,1), XX21(1,1), XX22(1,1), XX23(1,1), XX24(1,1)];
% %disp('    SkiprD=1.276      Unkpr=1.62       ShareSK=0.315  SharePRO=0.585   ShareSER=0.0991   Expshare=0.1  Patshare=0.61')
% Targets = {'SkiprD=1.276';'Unkpr=1.62';'ShareSK=0.315';'SharePRO=0.585';'ShareSER=0.0991';'Expshare=0.1';'Patshare=0.61'}
% disp(Moments)


% Targets: 1) SkillpreD=1.276, 2) UnskillpremD=1.62, 3) Share of skilled D= 0.315,
%          4) Share of workers in production = 0.58, 5) share of workers in services=0.091, 6) Patent Share=0.61

 
 Moments = [XX16, XX16b, XX18, XX20, XX21, XX24, XX22, betaPikD] %Benchmark: betaPikD is the US 1980 Corporate Equity/National Income, taken from Thomas Piketty & Gabriel Zucman "Capital is Back: Wealth-Income Ratios in Rich Countries 1700-2010" file USAxls. available at http://piketty.pse.ens.fr/fr/capitalisback
 DATA=[1.28, 1.37, 1.62, 0.315, 0.585, 0.63, 0.099, 0.67];
 

Endogenous = [the0D, the0F, wDL]

% penalties

p1 = 1;
p2 = 1;
p3 = 1;
p4 = 1;
p5 = 1;
p6 = 1;
p7 =1;
p8 =1;

pen_vec =[p1^2, p2^2, p3^2, p4^2, p5^2, p6^2, p7^2, p8^2];
Weight = diag((1./(DATA.^2)))*diag(pen_vec); % Weight matrix


Dist=(Moments-DATA);
Objective =(Dist*Weight*Dist') % Loss fuction: model-moments distance


disp('% matched')

 Fit = [XX16/1.276, XX16b/1.37 XX18/1.62, XX20/0.315, XX21/0.585, XX24/0.63, betaPikD/0.67,  XX22/0.099 ] % model fit


% XXstud
% Total_population


X